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Abstract. We perform global three-dimensional simulations 
of accretion disks integrating the compressible, non-viscous, 
but diffusive MHD equations. The disk is supposed to be 
isothermal. We make use of the ZEUS-3D code integrating 
the MHD equations and added magnetic diffusivity. We mea- 
sure the efficiency of the angular-momentum transport. Vari- 
ous model simulations delivered transport parameters of ass = 
0.01 to 0.05 which are consistent with several local numerical 
investigations. Two of the models reach a highly turbulent state 
at which ass is of the order of 0. 1 . After a certain stage of satu- 
rating of the turbulence, Reynolds stress is found to be negative 
(inward transport) in many of the models, whereas Maxwell 
stresses dominate and deliver a positive (outward) total trans- 
port. Several of the models yield strongly fluctuating Reynolds 
stresses, while Maxwell stresses are smooth and always trans- 
port outwards. Dynamo action is found in the accretion disk 
simulations. A positive dynamo-a is indicated in the northern 
hemisphere of the most prominent run, coming along with neg- 
ative kinetic and current helicities (all having the opposite sign 
on the southern side). The dipolar structure of the magnetic 
field is maintained throughout the simulations, although indi- 
cation for a decay of antisymmetry is found. The simulations 
covered relatively thick disks, and results of thin-disk dynamo 
models showing quadrupolar fields may not be compatible with 
the results presented here. 
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1. Introduction 

Accretion processes in astrophysical disks lead to enormous 
luminosity and huge changes in disk structure during astro- 
nomically short times. Efficient transport mechanisms are nec- 
essary to achieve such short time-scales. Anisotropic turbu- 
lence appears to be a major physical condition to provide astro- 
physical disks with strong transport. As these disks generally 
exhibit increasing specific angular momentum towards larger 
radii and thus fulfill the stability criterion of Rayleigh, they do 
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not give rise to an instability leading to turbulence by them- 
selves. Searches for instabilities in disks with rotation profiles 
similar to a Keplerian one unveiled several ways to turbulence 
being more or less favorable with respect to their prerequisites 
for the disk configuration. Gravitational instability needs the 
disk to be either cool or massive. Nonlinear and nonaxisym- 
metric perturbations require a severe additional perturber near 
the disk; conditions for instability in a purely hydrodynamical 
disk were derived by Dubrulle (1993). Hydrodynamic instabil- 
ity essentially comes down to violating the Rayleigh criterion 
saying that a rotation profile with an increasing specific angular 
momentum with radius is hydrodynamically stable, that is for 
dl 2 jdr > 0. Since the typical length of perturbations caused 
be external forces is supposed to be very large, the required am- 
plitude of the perturbations has to be considerable, too, in order 
to violate the Rayleigh criterion locally. If the length scale of 
the perturbation is comparable to the radius, a strong alteration 
of the Keplerian velocity of several per cent is needed. Con- 
vection was shown to deliver either negligible transport (Stone 
& Balbus 1996) or inward angular momentum transport (Kley 
et al. 1993). 

The requirements for the magnetic shear-flow instability 
(Balbus & Hawley 1991) do match astrophysical conditions in 
accretion disks in many configurations. All it needs is a radially 
decreasing angular velocity and a weak magnetic field thread- 
ing the rotating object. It can even be shown that the temper- 
ature range applicable to the magnetic shear-flow concept is 
very broad; even very small ionization fractions are sufficient 
to magnetize a disk in many cases (Balbus & Hawley 1998). 

First numerical approaches to the magnetic shear-flow in- 
stability tackled the local problem; the linear analysis as de- 
scribed by Balbus & Hawley (1991) were immediately fol- 
lowed by 2D simulations (Hawley & Balbus 1991) of a small 
box-shaped domain which was cut out of the disk. These 
computations confirmed the relation between magnetic-field 
strength and wavenumber derived from the linear analysis 
earlier, they showed that the maximum growth rate of any 
wavenumber is independent from the field strength and that 
the system is capable of transporting angular momentum. Be- 
cause of being restricted to axisymmetric problems, they could 
not provide self-sustained turbulence which needs dynamo- 
action, and the slow decay of the turbulence is an indirect effect 
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of the Cowling theorem. Improved computations dealt with a 
three-dimensional even though local configuration, and partic- 
ular care was taken for the radial boundary conditions which 
are not simply periodic, but account for the shear due to dif- 
ferential rotation (Hawley et al. 1995). These comprehensive 
computations indeed resulted in magnetically sustained turbu- 
lence whose anisotropy causes efficient outward angular mo- 
mentum transport. This work was followed up by stratified 
models (Stone et al. 1996) covering more than 50 orbital pe- 
riods of the box cut-out. As the computational domain covered 
2 density scale heights, this work was a first step towards global 
simulations, followed up by similar approaches such as Ziegler 
& Rudiger (2000a, b). 

Linear studies of global configurations of disks threaded by 
magnetic fields in various directions were carried out. Curry 
& Pudritz (1995) investigated the stability for vertical and az- 
imuthal fields threading the disk. They found in detail that the 
actual initial field geometry does not strongly depend on field 
topology as was suggested by the numerical simulations of 
Hawley & Balbus (1991). Rudiger et al. (1999) particularly ad- 
dressed the angular momentum transport in their linear study. 
These investigations are now followed by nonlinear simula- 
tions of a compressible fluid with density stratification in global 
computational domains. The most recent global approach was 
presented by Hawley (2000) following the evolution of a thick 
torus under the influence of an external magnetic field thread- 
ing parts of the computational domain. The magnetic shear- 
flow instability was found to set in quickly causing enough 
turbulent viscosity to soon form a Keplerian velocity profile. 
All the above mentioned numerical studies integrate the ideal 
MHD equations omitting magnetic diffusivity. Local simula- 
tions including diffusivity were performed by Fleming et al. 
(2000) which proved the onset of instability even for low con- 
ductivity. 

The full understanding of accretion disks implies self- 
excited dynamo action as well as angular momentum transport. 
Can positive angular momentum transport be coupled with a 
suitable kinetic helicity providing the expected dynamo action 
according to the a-effect principle? As the kinetic helicity in 
stratified, rotating disks is expected to be negative, a wrong 
sign (positive) would follow for the dynamo-a-effect accord- 
ing to the conventional a-theories. We present global three- 
dimensional diffusive simulations and study the angular mo- 
mentum transport as well as dynamo action and the sign of 
a dyn as a conse q Uence f the correlation with the flow. 

The present computations are not meant as particular sim- 
ulations of a star formation process. They focus on the appli- 
cability of the magnetic shear-flow instability for an astronom- 
ically fast process like the formation of a star. Such global sim- 
ulations are still a challenge for modern computers and fast al- 
gorithms. 

2. The simulation setup 

The computations presented here make use of the ZEUS-3D 
code developed for astrophysical problems of magnetohydro- 



dynamics (see the key papers Stone & Norman 1992a,b; Stone 
et al. 1992 for numerical, Clarke et al. 1994 for technical de- 
tails). We use cylindrical coordinates and an extensive compu- 
tational domain covering radii r = 5 to 6 and r = 4 to 6, 
a vertical extension of h = —1 to +1, and the full azimuthal 
range of cf> — to 2ir. See Table |l]for a list of models presented 
in this Paper. In this approach, we assume an isothermal disk to 
save computation time on the energy equation. The remaining 
system for integration is 

^+div(/m) = (1) 
d ptii 

+ div(puu) = — gradp — p grad $ + J x B (2) 
dB 

— = curlfw x B) + tjAB, (3) 
at 

where p, u, and B are the density, velocity, and magnetic field 
resp.; p is the pressure (p a p in our model), $ is the grav- 
itational potential (solely from a central mass M), J is the 
current density, and r\ is the magnetic diffusivity which is not 
an original ingredient to the ZEUS code. It is constant in time 
and space. The computations of the electromotive force in the 
routines emf s, mocemf s, and hsmoc of ZEUS-3D are ex- 
tended with the appropriate -qcxxAB components. We chose 
rj = 0.001 to 0.01. The additional time-step criterion resulting 
from the diffusivity is roughly 0. 1 for the finest grid used here. 
It is therefore irrelevant for the determination of the time-step 
which is typically 10~ 4 or, in the case of strong fluctuations 
of velocity and magnetic field, one or two orders of magni- 
tude lower. In physical units, this diffusivity is large and in fact 
accounting for a subgrid turbulent diffusivity which cannot be 
resolved. 

The gravitational potential is spherically symmetric, 
whence the z-component of the gravitation is retained within 
the disk. We therefore obtain a density stratification unlike the 
computations by Armitage (1998) who omits the z-component 
of the gravitation and applies periodic boundary conditions for 
the upper and lower boundaries. The sound speed is always 
c ac = 10 which is roughly 0.07uk in terms of the average 
Keplerian velocity in the simulated ring. The initial density 
distribution is Gaussian with r-dependent density scale-height, 
p(r,z) = p c exp(-z 2 /H(r) 2 ). 

The magnetic field threads the disk vertically with a mere 
z-component. We tried a homogeneous initial field B z (r) = 
const for the z-boundaries as well as a field 

B z {r, t = 0) = bor- 1 sin[27r(r - n )/(r - n )] (4) 

- where r; and r Q are the inner and outer boundary radii resp. - 
which has zero total magnetic flux through the z = const sur- 
faces. The choice of the initial field will have implications for 
the topology of the field later on. As the magnetic flux through 
the surfaces is kept constant, the first approach will always pre- 
serve a non-zero average magnetic field through the vertical 
boundaries, whereas the field penetrating the top and the bot- 
tom of the computational domain can completely vanish in the 
second choice of initial fields. The results presented here were 
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Table 1. Model parameter used for this study. Tq is the time when the magnetic field was switched on (in orbital units). The z- 
range is always from —1.0 to +1.0. If accretion boundary conditions occur, they refer to the outer boundary; the inner boundary 
is 'outflow' then. 



Model 


Resolution (z, r, (f>) 


i Radial boundary condition 




r -range 




To 


V 


la 


31 x 61 x 351 


inner outflow 




5.0-6.0 


0.222 


4.5 


0.001 


lb 


31 x 61 x 351 


all outflow 




5.0-6.0 


0.222 


4.5 


0.001 


11 


31 x 61 x 351 


accretion u- ln = — 0.001c ac 




4.0-6.0 


0.159 


9.4 


0.001 


III 


31 x 61 x 351 


accretion u- ln = — 0.01c ac 




4.0-6.0 


0.159 


41.2 


0.001 


IV 


31 x 31 x 351 


accretion ui n = — 0.001c ac 




5.0-6.0 


0.222 


26.8 


0.001 


V 


31 x 61 x 351 


Gaussian accretion U[ n — - 


0.001c ac 


4.0-6.0 


0.159 


29.6 


0.01 


VI 


31 x 61 x 351 


Gaussian accretion iti n = - 


0.01c ac 


4.0-6.0 


0.159 


18.9 


0.01 


VII 


31 x 31 x 151 


Gaussian accretion Ui n = - 


O.OOlCac 


4.0-6.0 


0.159 


16.1 


0.01 


VIII 


31 x 61 x 351 


accretion u- ln = — 0.001c ac 




3.0-7.0 


0.103 


24.7 


0.001 



obtained with the second approach of vanishing average initial 
field B z . The parameter b was chosen between 50 and 100 giv- 
ing an amplitude of the initial magnetic field of max(i3 z ) = 10 
to 20 corresponding to maximum Alfven speeds at the upper 
and lower boundaries of ua = 13 to 26 or ua = 1.3 c ac to 
2.6 c ac . The Alfven velocity in the equatorial plane is subther- 
mal with u A = 0.32 to 0.63 or 0.032 c ac to 0.063 c ac . 

The initial velocity field is a merely Keplerian motion fol- 
lowing = yj GM/r, where G is the gravitational constant 
and M is the central mass which is 10 5 in our computational 
units. The total disk mass is 35,000 for Model la, lb, and IV 
and 56,000 for the other models. Times are henceforth mea- 
sured in orbital periods which convert by T or b = 0.159 from 
the arbitrary units of the code. 

The number of cells in each coordinate direction was 31 x 
61 x 351 for the z-, r-, and (^-directions. Tests with up to 621 
cells in (/(-direction have been carried out, but the increased 
computation time does not allow reasonable periods to be cov- 
ered by the simulation. The aspect ratio of the grid cells is not 
unity. Finest resolution is achieved in radial direction, lowest 
in azimuthal direction: Az/ Ar = 2, rA<j)/Ar = 2.18 for the 
smallest radius and rA(j>/Ar — 3.28 for the outer radial edge. 
The fastest growing mode of the magneto-rotational instabil- 
ity has a wavelength of Ai ns t = 2^-^/16/15 ua/^- With the 
angular velocity of Cl = 28.2 in the middle of the computa- 
tion domain, we obtain wavelengths between 0.4 and 0.9 with 
the above given Alfven velocities of the initial field strength. 
These wavelengths are upper limits as they refer to the sites of 
maximum field and minimum density which need not coincide 
necessarily. Yet, they are well resolved by the computation do- 
main covering z -ranges of 2.0 and two different r-intervals of 
1.0 and 2.0. 

The integration of the magnetic fields was done with the 
evolution of the electromotive forces by the Constraint Trans- 
port scheme (CT, Evans & Hawley 1988) which preserves a 
divergence-free configuration. 

Two terms adding a source of viscosity to the hydrody- 
namic equations extend the above Navier-Stokes equation. The 
von Neumann-Richtmyer artificial viscosity depends on the 
square of velocity gradients and acts only on decreasing ve- 
locity in the direction of propagation, i.e. compression. The 



strength of this term is denoted by q con in the implementation 
of the ZEUS-3D code. The second term depends linearly on the 
velocity gradient and acts on both compression and expansion; 
the strength is represented by qn n . Tests with a stable Keplerian 
flow found the choice of g con = 0.1 and q\ in = 1.0 to be suit- 
able. These values were used throughout all the runs described 
here. The Courant number determining the 'safety' of a certain 
maximum allowed time step derived from the velocities, mag- 
netic fields, and the artificial viscosities, is set to 0.5. 

The ZEUS-3D code provides three methods of interpola- 
tion in the transport step: the first-order donor-cell method, the 
second-order van-Leer method, and the third-order piecewise- 
parabolic-advection method. The original ZEUS-3D copy 
chooses the second-order interpolation for the transport of all 
quantities, and we have not altered this switch for the compu- 
tations given here. 

Occasionally, pockets of extremely low density may 
emerge coupled with extraordinarily high speeds exceeding the 
Keplerian velocity. We suppress the existence of such pockets 
by a mass replenishment as soon as the density of a certain cell 
drops below 10~ 4 the central density. The actual mass being 
thus added to the total disk mass is found to be negligible. The 
mass replenishment thus acts like a lower limit for the time 
step. 

2.1. Boundary conditions 

The vertical boundary condition on the z = const-faces of our 
computational domain are closed for the flow, i.e. v z = 0. The 
boundary condition is stress-free in the sense of dv r /dz = 
and dv^/dz = 0. The magnetic field has to fulfill the 'op- 
posite' boundary condition as it is only allowed to cross the 
boundary perpendicularly, which is actually implemented af- 
fecting the electromotive force £ = u x B — r\ curl B such 
that E z = 0, d£ r /dz = 0, and dS^/dz = 0. 

No magnetic field lines are allowed to penetrate the radial 
boundary surfaces. Computations were carried out with two 
choices of radial boundary conditions regarding the velocity. 
The first does not allow flow through the top, bottom and outer 
radial boundaries. The inner boundary was chosen to be 'out- 
flow', that is, the physical arrays are constantly extrapolated 
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into the ghost zones of the computational domain. The only 
exception is an inflow, in our case a u r > at the inner bound- 
ary. This velocity component is reset to zero then. All radial 
velocities u r < are copied to the ghost zones beyond the 
computational domain. We will henceforth refer to these mod- 
els as Model la (only inner boundary has 'outflow') and lb (all 
both radial and both vertical boundaries have 'outflow'). 

Since the outflow condition is likely to empty the disk 
on the viscous time-scale, the second choice tries to account 
for the accreted matter and feeds the disk at the outer radial 
boundary. We sum up the mass loss rate at the inner bound- 
ary M = r min A(j>Az J2i J2k Pik{u r )ik for the total mass loss. 
A constant, small inflow velocity u- m is assumed for the outer 
boundary. Accordingly, the density at the outer boundary is de- 
termined to account for the mass flux at the inner boundary. 
Since we use the average influx, its z-dependence (nor even 
the 0-dependence) is not used for a z- (nor even </>-) depen- 
dent inflow from far radii. We tested inflow conditions with ho- 
mogeneous and Gaussian and density distribution. It was sup- 
posed that an average influx of mass will establish according to 
the current viscosity (either numerical or turbulent) in the disk. 
These models are refered to as Model II- VIII in the following. 

The azimuthal velocity, u^, is extrapolated into the radial 
boundary zones by a power law r -15 based on the last zone of 
the computational domain. The velocities are thus not Keple- 
rian, they just follow the same radial power law. Since we cover 
the full azimuthal range in <j>, the boundaries at the <f> = const 
surfaces are periodic, naturally. 

Even though the ZEUS-3D code is remarkably stable for a 
broad variety of configurations and physical problems, a com- 
prehensive study of the influence of the computational schemes 
provided were found to be inevitable for trustworthy results. 

3. Patterns, transport, spectra 

3.1. Glancing at the solutions 

Typical vertical slices from a run of Model II are shown in 
Figures [j] and |[ The central mass is located on the left side, 
outside the computational domain. The shadings represent the 
azimuthal components of the magnetic and velocity fields, 
while the arrows are the (r, z) -projections of the actual vectors. 
Shaded areas near the equatorial plane of the disk show the typ- 
ical Keplerian velocity profile. However, at high altitudes above 
the equator, a strong super-rotation emerges coupled with rela- 
tively high infall velocities. (These high radial velocities have 
been cut in Figure [l] to better show the velocities in the disk.) 

The computational domain thus covers part of the corona of 
the disk where the density falls below 1CP 3 times the equatorial 
density. The infall velocity in the corona is of the same order 
of magnitude as the Keplerian velocity, whereas the non-orbital 
velocities within the disk are much smaller. The actual flow in 
the disk as shown in Figure |l] is limited by 0.35c ac . 

After saturation of the turbulence, the temporal behavior of 
ass changes. The fairly smooth function turns into an oscilla- 
tory behavior. The typical feature of this moment is the zero- 
crossing of the Reynolds stress which is essentially negative all 
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Fig. 1. Vertical cut through the disk of Model II after 15.7 or- 
bital revolutions, excluding the high-velocity coronal compo- 
nent of the disk. The shading represents the azimuthal com- 
ponent, the arrows are the velocity vectors projected onto the 
(r, z)-plane 
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Fig. 2. Vertical cut through the disk of Model II after 15.7 or- 
bital revolutions with an azimuthal-component shading and 
poloidal magnetic-field vectors 



through the saturation parts of the simulations. It is expected 
that the infall of matter at the outer radial boundary tends to per- 
form relaxation motions, in particular for Models II-IV where 
matter falls in homogeneously. The exclusion of the outermost 
six cell planes from the computation of the spatially average 
ass does not alter the result (the actual relaxation motions are 
observed only in the outermost three cell planes). Relaxation 
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motions will consist mainly of u z components which do not 
affect radial angular-momentum transport directly. 

A problem connected with the homogeneous feeding of the 
disk from the outer boundary is the appearance of a 'density 
front' . This is the reason why the density at the outer boundary 
is highest which may look odd in the context of an accretion 
disk. This incoming density enhancement can be partly leveled 
by the presumed infall velocity at the outer boundary. 

3.2. Angular momentum transport 

According to the standard model of a turbulent disk accord- 
ing to Shakura & Sunyaev (1973), the r0-element of the stress 
tensor of velocity and magnetic-field fluctuations, 



W riP = 



(P) 



scales like the speed of sound such as 

W r A> = a ss cl 



.2 

-ac ' 



(5) 



(6) 



parameterized by an unknown agg which is smaller than unity 
in the case of subsonic turbulence. We normalize the stress by 
the average density in the computational domain; other authors 
used the density averaged in the equatorial plane leading to 
lower limits for agg. This stress measures the angular momen- 
tum transport in the disk; positive values mean outward trans- 
port. We average this quantity over the entire computational 
domain according to 



/ pv! r u'^rd4>dzdr 
J rd(j)dzdr 

for the velocity and 

1 / B'rB'^rdfidzdr 
Ait J rd(j)dzdr 



(7) 



(8) 



for the magnetic field fluctuations. As the azimuthal velocity 
contains a large supersonic axisymmetric mode - the Keplerian 
flow - we adopt u'^ = u^,— ^/GMjr for the azimuthal velocity 
fluctuations. The averaging is in fact carried out on a discrete 
computational mesh by 



ink E 



i'3 



for the velocity part and 

1 E,-^Ei T,k B 'r B 'd> 



4tt 



n %nk Ej r 3 



(10) 



for the magnetic part. 

The time series of the outflow model la is shown in Fig. ^|. 
Model la in particular shows the typical short-lived occurrence 
of an equivalent to the channel solution as was found in lo- 
cal box simulations; a similar but less pronounced behavior is 
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Fig. 3. Angular momentum transport by Reynolds (dashed) and 
Maxwell (dotted) stresses as derived from a run of Model la, 
that is with outflow boundary condition at the inner radial face 
and a closed boundary for flow and magnetic fields at the outer 
radial face. The total transport (thick solid line) is positive and 
of the order of 0.0005. 
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(9) pjg_ 4_ Energy in individual velocity components of Model la 



found in Model lb. This period is the only time in both Mod- 
els la and lb when Reynolds stresses dominate the Maxwell 
stresses. 

Time series for long runs of Models II and Model III 
are shown in Figure |5] and [7j Generally, the accretion model 
shows stronger angular momentum transport than the the sim- 
pler inner-outflow Model la. Models la and lb run out of matter 
after about 8 orbital revolutions with enabled magnetic fields, 
and the agg values are no longer meaningful. 
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Fig. 5. Angular momentum transport by Reynolds (dashed) and 
Maxwell (dotted) stresses as derived from a run of Model II. 
The total transport (thick solid line) is positive and of the order 
of 0.01 



Fig. 7. Angular momentum transport by Reynolds (dashed) and 
Maxwell (dotted) stresses as derived from a run of Model III. 
The homogeneous inflow velocity is ten times higher in this 
simulation that in Figure || 
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Fig. 6. Energy in individual velocity components of Model II 



A striking feature of the accretion Models II and III is the 
change of sign in the Reynolds stress after about 4-6 orbits. 
We find that, in a long run, the outward angular momentum 
transport is a sole consequence of the magnetic stresses. Addi- 
tionally, the Reynolds part shows strong amplitude variations 
unlike the Maxwell stress. 

Figures [| and ^ show the respective temporal evolution of 
the fluctuative kinetic energies. For the z- and r-components, 
the zero mode is subtracted such as (v! 2 ) = {u 2 z ) — (u z ) 2 \ 
the azimuthal fluctuative energy is found by (u'£) = {{u^ — 



y^GM/r) 2 ). All three models show a clear energy 'sorting' of 
the z-, r-, and (^-components in this order. 

A particularly long run was performed with a computa- 
tional domain only half as large as Models II and III and is 
denoted by IV. The initial and boundary conditions are like 
in Model II. The long-term behavior is not too similar to that 
of Model II with very strong oscillations of Reynolds versus 
Maxwell stresses. A behavior comparable with Model II may 
be spotted during the first 8 orbital periods. After that time, the 
energy sorting of the radial and vertical components is given 
up. All models from la to IV appear to be not fully saturated 
with respect to the ever-growing kinetic energies. 

We should note that a homogeneous infall might be suitable 
for a realistic disk in the sense that the infalling matter does 
not 'know' the structure of the disk. Nevertheless, numerical 
problems may occur when the incoming homogeneous 'matter 
front' meets the near-Gaussian density structure. Models V and 
VI apply an inflow with pre-defined Gaussian infall profile over 
z. The respective ass-profiles are shown in Figures |and[l0[ 
the first having an infall velocity of u la = —0.01, the second 
""in = —0.1. The initial behavior is similar to the respective 
Models II and III including the change of sign of the Reynolds 
stress after a couple of orbits. After as many as 10 orbital pe- 
riods, the disk turns into a significantly different regime with 
enhanced outward transport in both models. 

The respective kinetic energies are given in Figures ^ and 
pTj In contrast to the above studied models, energies do not de- 
velop in an ever-growing way, but appear to fluctuate about 
an average. The only exception is the regime change near 
t = 42 T orb (Model V) and t = 29 T orb (Model VI) when 



R. Arlt & G. Rudiger: Global simulations of accretion disks 



7 



the average (essentially of (uf)) changes, but the stationary 
appearance remains. 
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Fig. 8. Angular momentum transport by Reynolds (dashed) and 
Maxwell (dotted) stresses as derived from a run of Model V. 
Density is not falling in homogeneously, but with a predefined 
Gaussian vertical distribution; the (homogeneous) inflow ve- 
locity is 0.01 
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Fig. 9. Energy in individual velocity components of Model V 

The distribution of ass versus time and the z-direction is 
shown in Figure [l2| for Models V. Outward angular momen- 
tum transport is marked with white areas, inward transport by 
dark areas. Strong transport is observed at large distances from 
the equator shortly after the onset of the instability. Regions of 
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Fig. 10. Angular momentum transport by Reynolds (dashed) 
and Maxwell (dotted) stresses as derived from a run of Model 
VI. The inflow velocity is ten times higher than in model V, 
Figure || 
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Fig. 11. Energy in individual velocity components of Model VI 



outward transport migrate towards the equator. This is particu- 
larly marked in the Model-V plot for which we found a devel- 
opment into a strong-transport regime (cf. Figure |]). Positive 
angular momentum transport is not restricted to high altitudes 
anymore. 

3.3. Spectral decomposition 

Power spectra resulting from Fourier transforms of the mag- 
netic field of Model V over the the azimuthal and radial di- 
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Table 2. Results of the models given in Table 111. The duration is 
the total period of the run with enabled magnetic fields. The last 
column gives the period from which the ass" avera g e was ob- 
tained. Times are in units orbital revolutions of the inner edge 



TIME = 37.75 orbits 



Mod. 


Resolution 


Duration 


(ass) 


during 


Ia 


31 x 61 x 351 


8.3 


0.0003 


8.5-13.4 


lb 


31 x 61 x 351 


12.2 


0.01-0.02 


8.0-13.5 


II 


31 x 61 x 351 


14.7 


0.014 


13.0-24.0 


III 


31 x 61 x 351 


9.7 


0.041 


44.0-50.9 


IV 


31 x 31 x 351 


41.5 


0.014 


35.0-69.8 


V 


31 x 61 x 351 


11.9 


0.038 


33.0-41.5 








0.12 


42.0-45.5 


VI 


31 x 61 x 351 


16.1 


0.037 


22.0-29.0 








0.14 


30.0-35.0 


VII 


31 x 31 x 151 


36.8 


0.080 


32.0-44.0 


VIII 


31 x 61 x 351 


22.4 
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Fig. 12. Distribution of the total stress ass versus time and 
vertical coordinate for Model V. Light areas represent inward 
transport, dark areas outward transport 

rections are shown in Figures and [l4], respectively. The az- 
imuthal decomposition of the velocity field is shown in Fig. [l5|. 
The individual azimuthal spectra for all points in the (z,r)- 
plane were averaged, as were the points of the ((f), z)-plane for 
the radial spectra. In the azimuthal spectra, we observe slight 
maxima in power at k = 5 and k = 8, followed by a strong, 
roughly power-law decline between k ~ 20 and k ~ 100. The 
magnetic spectrum is steeper than the kinetic one. The under- 
lying power-laws would be E oc fc~ 3 to E oc k~ A which is 
significantly steeper than a Kolmogorov spectrum, and is sim- 
ilar to what Armitage (1998) found for a global simulation in 
the far-fc range. The Kolmogorov spectrum of isotropic turbu- 
lence would exhibit a wavenumber exponent of — 5/ 3 ; here we 
face MHD turbulence which implies an effect of Alfven waves 
impeding the transfer of energy towards smaller scales. The 
energy spectrum of isotropic MHD turbulence is E oc fc -3 / 2 , 
whence more shallow than the Kolmogorov spectrum. Con- 
trasting with the azimuthal decomposition, the radial spectrum 
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Fig. 13. Azimuthal Fourier decomposition of the velocity com- 
ponents after 37.75 orbital periods in Model V 
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Fig. 14. Radial Fourier decomposition of the magnetic-field 
components after 37.75 orbital periods in Model V 



matches the theoretical MHD spectrum over the entire range of 
modes as can be seen in Fig. [14]. 

3.4. Accretion rate 

Choosing the rough dimensions of an inner protostellar system 
of 100 AU radius for scaling the mass flux through our mod- 
els, we obtain an accretion rate in physical units as shown in 



Fig. 16 Widely accepted values for the accretion rates of pro- 
tostellar disks around solar-type stars are 10~ 7 to 10 -6 M Q /yr. 
The compilation of results for those stars by Strom (1994) 
rather favors 10 -6 M Q /yr. The accretion rate increases with 
Stella mass and may already exceed 10~ 5 Af Q /yr for a 3-solar- 
mass star. Our results such as given in Fig. |l^ are thus not too 
far from the accretion rates deduced from the observations. 



R. Arlt & G. Rudiger: Global simulations of accretion disks 



9 



TIME = 37.75 orbits 



10 L 



10 



« 

S= 10 4 
o 



1CT 




u y (m=0: 8.5E+10) 
u r (m=0: 2.0E+07) 
u, (m=0: 3.8E+06) 

10 

AZIMUTHAL MODE 



100 



Fig. 15. Azimuthal Fourier decomposition of the velocity com- 
ponents after 37.75 orbital periods in Model V 
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Fig. 16. Accretion rate of Model II scaled to a system at 1 00 AU 
distance from the central object. 



3.5. Comparison with other studies 

Among stratified local simulations, we find a number of very 
similar results as regards the magnitude of ass, sucn as m 
Stone et al. (1996) who obtained 0.01, Hawley et al. (1996) 
where an average of 0.016 was found, with typical values of 
0.01. The largest box used delivered 0.02, the smallest, 0.003. 
Upon studying the dependence of the turbulence on the rotation 
law fi oc r~ q , Hawley et al. (1999) found 10" 2 for q = 1.5. A 
most recent local-box simulation by Ziegler & Rudiger (2000b) 
covering 288 orbital periods resulted in a space-time average of 
0.015, where the Reynolds part was 2.8 • 10~ 3 and the Maxwell 
part was 1.2 • 10~ 2 . The stress was normalized with the aver- 
aged equatorial-plane pressure and thus provides a lower limit 
for ass- The early stratified-box simulations by Brandenburg et 
al. (1995) provided ass an order of magnitude lower between 
0.001 and 0.005. All the models agree upon the dominance of 
Maxwell stresses over Reynolds stresses. 



It is most interesting to note that other global simulation 
such as performed by Armitage (1998) or Hawley (2000) re- 
sult in an order of magnitude higher values of 0.2 to 0.3. Con- 
trasting with this, the global spherical, but unstratified simula- 
tion by Drecker (2000) gives again an ass which is an order 
of magnitude lower than shown here, i.e. 10~ 3 . Those simu- 
lations alone included a physical viscosity. The simulations of 
Hawley (2000) - and with somewhat different objectives, those 
of Machida et al. (2000) - had open boundaries in radial and 
vertical directions. This fact limits the run-time of the models 
which are prone to running out of matter. 



4. Dynamo effect 

The long-term existence of non-zero toroidal fields in each 
hemisphere as shown for example in Figure |l9| for Model V 
indicates the possible action of a dynamo in the disk. 

The local -box simulations by Hawley et al. (1996) impose 
severe constraints on the growth of an average magnetic field 
because of the periodicity of the boundary conditions. A non- 
zero term of the average magnetic-field growth is due to the 
shear and requires the presence of an average radial component, 
in the local nomenclature (B x ). We follow their derivation of 
the averages by replacing the volume curl integral by a vector- 
produc surface integral of (3). The average 



d(B) 
dt 



- \ - / cml(u X B - r>cvx\B)dV (11) 
?rAr 2 Az J y 1 ' 



is equivalent to 

d(B) 
dt 



T-=-r / dS X (ux B -t? curl B) (12) 
TrAr z Az J 



and delivers average components as 



d(B z ) 




dt 


\ dr 


d(B r ) 


(~ 


dt 


\dz 


d(B^) 


(- 


dt 


\dz 




(- 




\ dr 



u r B z ) + 77 

II,, B; ) - 1] 

V 



d dB z 
dz dr 
1 d dB z d 2 B d 



r dz dcj) 
d ldrBj 



dz 2 



dr r dr 



(13) 
(14) 

(15) 



These considerations show that slight initial fluctuations in B 
are able to create non-zero average magnetic fields at any time, 
even if their (B 1 ) vanishes, as a consequence of shear and ra- 
dial accretion through u ul . In particular, the non-constant ver- 
tical profile of the orbital velocity produces a B^ from the ini- 
tial perturbation B z . Fluctuations in the accretion flow generate 
B r and B z . The actual radial differential rotation then quickly 
produces B^ from even very small B r , and toroidal fields are 
prone to dominate the disk fields. 
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4.1. Sketching mean-field dynamos 

Kinematic dynamos have been extensively examined in numer- 
ous previous publications concerning various types of geome- 
tries. A widely used mean-field approach to circumnavigate the 
difficulties with resolving the small scales in a model applies 
the likely amplification of a magnetic field through a fluctuative 
electromotive force parallel to the actual large-scale magnetic 
field. This effect is generally expressed by the term 'a-effect' 
and can be thoroughly studied in e.g. Krause & Radler (1981). 
In a model splitting large and small scales, the induction equa- 
tion extends beyond the large-scale electromotive force such 
as 



d(B) 
dt 



= curl«u) x (B)+S). 



(16) 



The a-effect is part of the development of the small-scale elec- 
tromotive force as 



Si 



dyn 

a ■ 

dyn 



(Bj) - rjijk(B k ) tj + .... 



(17) 



where a t J and rjijk are tensors in general. If written in compo- 
nents, the induction equation shows that the essential compo- 
nent of the a dyn tensor is the a^ n component which converts 
azimuthal fields into radial and vertical fields, whereas the dif- 
ferential rotation converts radial and vertical fields back into 
the toroidal component by the large-scale part of (u) x (B). 
This approximation is usually referred to as an ail-dynamo. 

Second-order correlation approximation leads to a relation 
of a dyn with the helicity of the flow, 



oo 

a dyn = -\J { u \t) ■ curlu'(t - t)} dr, 



(18) 



which is, when approximated by a typical time-scale r corr , 



v d y n 



(u' ■ curl it') Tr, 



(19) 



The same principle is applicable to the current helicity (Keinigs 
1983) giving 



, d yn 



-^(B'-curlB')- 



(20) 
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Fig. 17. Total a ss and parity for Model II. A parity of -1 de- 
notes fully antisymmetric (dipolar) solutions 
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Among other issues, the evolution of these quantities for the 
helicity will be subject of the following sections. 



Fig. 18. Total a S s and parity for Model V. 



4.2. Symmetry of magnetic fields 

Special attention is payed to the symmetry of the magnetic 
fields of the solutions. In general, the preference of a certain 
symmetry is a tool to check the consistency with mean-field dy- 
namo models. Additionally, dipolar (antisymmetric) fields are 
supposed to support the launch of winds forming jets from the 
disk better than quadrupolar (symmetric) fields, with the latter 
rather providing closed field lines within the disk. The mean- 
field disk dynamo of Rekowski et al. (2000) produced dipolar 
magnetic fields for negative a dyn in the northern hemisphere. 

A number of mean-field dynamo models in one to three di- 
mensions delivered varying symmetry of the solutions with a 



major influence of the boundary conditions - vacuum versus 
perfectly conducting. Three-dimensional investigations of sta- 
bility by Meinel et al. (1990) and Elstner et al. (1992) found 
quadrupolar solutions for the condition of low-conductivity 
surroundings of the disk, whereas perfectly conducting bound- 
aries delivered dipolar solutions. A local model with periodic 
boundaries will thus not be able to answer the question of the 
final symmetry of the solutions. An ideal global model would 
comprise the entire investigated object and is less depending on 
the surroundings (which will be vacuum in very good approxi- 
mation). 
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The parity of a the toroidal magnetic field is measured by 

P= E , (2D 

where E s and i? a are the symmetric and antisymmetric energy 
parts resp., and E is the total energy of in toroidal field com- 
ponent. Figures |l7| and [l8| show the temporal development of 
parity of the toroidal magnetic fields in the Models II and V 
along with a repeated plot of ass- 

We should emphasize that the initial perturbation of the 
magnetic field had dipolar symmetry. Since the total flux 
through the vertical surfaces, however, vanishes, field config- 
urations with either symmetry may arise from the evolution of 
the model. The dipolar initial perturbation may not have been 
fully re-organized as the diffusive time-scale 

Tdiff = H 2 /rj (22) 

is 1600 rotation periods for the low-diffusivity model II and 
160 revolutions for the high-diffusivity model V. The latter 
covered 0.12 r<iig. An additional experiment used a completely 
mixed initial magnetic-field perturbation according to a distri- 
bution which has zero parity in B z , and the toroidal component 
immediately emerging from the slight vertical has zero parity, 
too. Again, the magnetic flux through the vertical surfaces is 
zero. Such initial conditions have been applied to Models II 
and V. The runs lasted roughly 10 revolution resulting in parity 
variations with a slight tendency towards quadrupolar (sym- 
metric) solution (P from +0.1 to +0.2). 

4.3. The dynamo of the disk model 

Average kinetic helicities were computed based on a layer 
slightly below one density scale-height in both hemispheres. 
The average is taken in <f>- and r-direction. There is a consider- 
able scatter in the time series of the helicity, but temporal av- 
erages are all negative on the northern side for all models. The 
only exception is the northern part of Model II with a near-zero 
value. Negative sign also holds for the two test models with 
mixed-parity initial perturbation. Average current helicities are 
also negative on the northern side (and positive on the south- 
ern one) throughout almost all models. The only exceptions 
are the mixed-parity models which show a current helicity with 
opposite sign compared with the kinetic helicity in both hemi- 
spheres. 

Since the angular momentum transport is dominated by 
magnetic stresses, we also ask about the energy in the flow 
fluctuations and the magnetic field. The temporal evolutions of 
these two quantities in Figures ^ and |T] show a clear differ- 
ence between the Models II and V. Kinetic-energy dominance 
holds for Model II, whereas magnetic dominance is found for 
Model V. 

In order to evaluate the applicability of the a dyn approach, 
we compute the average toroidal field (B^) and compare it with 
the actual average electromotive force 8^ — ((u) x {B))$ 
derived directly from the simulated vector fields. A correla- 
tion between the two quantities may justify the a-effect for 



dynamo generation of magnetic fields. As a meaningful a dyn 
must change its sign at the equator, we plot the correlation 
for both hemispheres of the disk separately; the result is given 
in Figure ^3|. As the temporal fluctuations are strong, we plot 
the time-averages of B^, and (u) x (B)^ in Figure ^2| only 
with their resulting sign. The quantities are also averaged in 
azimuthal direction and plotted in the (r, z)-plane. White ar- 
eas denote a negative sign; the dominance of negative averaged 
electromotive forces indicates a positive a dyn in the northern 
hemisphere in accordance with Figure 

The classical understanding of the Coriolis force giving 
preference to left-handed helicity on the northern hemisphere 
(whence positive a dyn ) appears not applicable for the corre- 
lation plots of Models II and III. The fluctuations were strong, 
and no meaningful sign of a dyn can be derived. However, Mod- 
els la, V, and VI show a mostly negative averaged EMF for 
both hemispheres, while (Bs) changes its sign. An indication 
for positive a dyn is thus found from these simulations. 

In fact, Models V and VI are those which show saturated ki- 
netic energies. The other significant difference to the indefinite 
Models II and III is the lOtimes higher diffusivity, r\ = 0.01. 
Notwithstanding, an indication for a negative a dyn was already 
found in simulations of a local box cut out of the disk by Bran- 
denburg et al. (1995). The same sign is found by Ziegler & 
Rudiger (2000b) from long runs of local simulations, although 
a clear influence of resistivity on the evolution was found 
there. Lowest magnetic Reynolds numbers are Re m = 1450 
for Model V in this Paper, calculated with the velocity differ- 
ence between inner and outer radial boundary, just as Ziegler 
& Rudiger did for the local box. 

These considerations are somewhat limited, since a 
straight-forward regression line in Fig. |23| will have a signif- 
icant offset from the origin of the graph; a vanishing (B^) does 
apparently not coincide with a vanishing ((it) x (B))^. The ef- 
fect of the other field components is considered small though; 
while the average hemispherical (B^) are of the order of 100, 
the (B r ) is of order 10 and (B z ) of order unity. A consider- 
able contribution is suspected from the current density through 
the T^jfc-term in (fl7|), and is appears to be quite natural that an 
offset in the correlation graphs in Figure ^3] is found. For this 
reason, we did not plot such regression lines in the Figure. 

Additionally, the sole consideration of the induction equa- 
tion in the "classical" dynamo theory leads to the omission of 
the Lorentz force, which is not only essential for the onset of 
the magneto-rotational instability but for the maintenance of 
the turbulence and the transport of angular momentum as well. 
It is thus obvious that the neglect of the back-reaction of mag- 
netic fields on the flow need not lead to representative models 
of Keplerian disks. 

5. Conclusions 

Three-dimensional global simulations have shown, that the 
magnetic shear-flow instability is a fast mechanism to gener- 
ate a turbulent flow in a Keplerian disk. We presented models 
with a computational domain covering the full azimuthal range, 
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Fig. 19. Average toroidal magnetic field of Model V. The 
dashed line is the toroidal field averaged only in the northern 
hemisphere, the dotted line is from the southern hemisphere 



Fig. 21. Temporal evolution of kinetic (solid line) and mag- 
netic (dashed line) energies of Model V. The energy is mea- 
sured as a radial and azimuthal total in a horizontal layer at 
z = +0.39 which is slightly below one density scale-height 
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Fig. 20. Temporal evolution of kinetic (solid line) and mag- 
netic (dashed line) energies of Model II. The energy is mea- 
sured as a radial and azimuthal total in a horizontal layer at 
z = +0.39 which is slightly below one density scale-height 



accounting for accretion of matter. The Shakura-Sunyaev pa- 
rameter ass increases rapidly during the first revolutions af- 
ter switching on the magnetic field and reaches 10~ 2 -10 _1 at 
this stage of the computations, though it appears not yet fully 
saturated in several of the models. Maxwell stresses exceed 
Reynolds stresses almost entirely. The latter undergo strong 
fluctuations and are partly negative. 



Indications for a dynamo action in the disk are found, the 
corresponding dynamo-a dyn tending to be negative north of the 
equatorial plane and positive south of the equatorial plane. The 
generated magnetic fields may maintain the turbulence even if 
the external field ceases. The strong excitation of low-order az- 
imuthal modes in the magnetic field is another promising fact 
for dynamo action with respect to the Cowling theorem. 

5.1. The dynamo of Model V 

For various reasons, a representative example for an accretion 
disk dynamo is provided by Model V (and nearly as well by the 
similar Model VI with higher inflow velocity). The time series 
is sufficiently long covering more than 18 revolution periods. 
The angular momentum transport reaches a saturated level of 
high efficiency during the last four orbital periods. The mag- 
netic diffusivity is roughly two orders of magnitude higher than 
the numerical diffusivity. The striking dominance of one sign of 
the averaged EMF in Figure |2| indicates a physically evolved, 
relevant state of the system. 

The results connected with the outcomes of mean-field dy- 
namo theory include include the following facts: (i) a negative 
average toroidal magnetic field is found for the northern hemi- 
sphere, a positive on the southern one. The averaged EMF is 
negative in both hemisphere indicating a positive a dyn -effect 
on the northern side (negative on the southern side), (ii) The 
correlation plot of average toroidal field and EMF gives the 
same picture, (iii) negative kinetic and magnetic helicities on 
the northern hemisphere (positive on the southern one) are also 
consistent with a positive a dyn . (iv) The dipolar structure of 
the solution contrasts with the results from mean-field dynamo 
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Fig. 22. Time-average toroidal magnetic field (TOP) and the 
toroidal electromotive force (BOTTOM) shown in spatial dis- 



models of disks which yield quadrupolar solutions for positive 
a Ayn . We have to add though that the disks of our Paper are 
not thin. Mean-field simulations by Covas et al. (1999) show a 
transition from quadrupolar to dipolar fields when the opening 
angle of the disk is enlarged, rather representing a torus than a 
disk. 

The extension of the run-times of such simulations exceed- 
ing the order of diffusion times promises further results about 
dynamo action in accretion disks. The connection with mean- 
field concepts are most interesting as well as the implications 
for jet-launch models. 
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Fig. 23. Correlation of the average electromotive force versus 
the mean B^, averaged for the northern and southern hemi- 
sphere separately. The quantities were derived, from top to bot- 
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correlation gives an indication for the sign of the dynamo a- 
effect. 
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